# Replication data set: Landholding Inequality, Social Control, and Mass Opposition to Suffrage Extension

library(sf)
library(AER)
library(lfe)
library(sandwich)
library(lmtest)
library(ivmodel)
library(texreg)
library(tidyverse)
library(modelsummary)
library(GGally)

##### Data

# Pooled Data

load("C:/Dropbox/Papers/Suffrage_Vote/Paper/Submission_BJPS_Final/pooled.RData")

# Complete (unpooled) Data

load("C:/Dropbox/Papers/Suffrage_Vote/Paper/Submission_BJPS_Final/all.RData")

# Different versions of the data with missings for some variables

load("C:/Dropbox/Papers/Suffrage_Vote/Paper/Submission_BJPS_Final/res1.RData")

load("C:/Dropbox/Papers/Suffrage_Vote/Paper/Submission_BJPS_Final/res2.RData")

load("C:/Dropbox/Papers/Suffrage_Vote/Paper/Submission_BJPS_Final/res3.RData")


############# Main Paper

##### Table 1: Landholding Inequality and Support for Suffrage Extension

summary(m1 <- felm(yessh ~ bzineq| vote | 0 | BEZIRK, data = pooled))

summary(m2 <- felm(yessh ~ cath+ger+urban+log(pop)+bzineq+splag| vote | 0 | BEZIRK, data = pooled))

# First stage nutrient levels  

summary(m31 <- felm(bzineq ~ nut| 0 | 0 | BEZIRK, data = pooled))

summary(m32 <- felm(bzineq ~ cath+ger+urban+log(pop)+nut+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage nutrient levels 

summary(m4 <-felm(yessh ~ 1| vote | (bzineq~nut) | BEZIRK, data = pooled))

summary(m5 <- felm(yessh ~ cath+ger+urban+log(pop)+splag| vote | (bzineq~nut) | BEZIRK, data = pooled))

# First stage wheat

summary(m61 <- felm(bzineq ~ whe| 0 | 0 | BEZIRK, data = pooled))

summary(m62 <- felm(bzineq ~ cath+ger+urban+log(pop)+whe+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage wheat

summary(m7 <- felm(yessh ~ 1| vote | (bzineq~whe) | BEZIRK, data = pooled))

summary(m8 <- felm(yessh ~ cath+ger+urban+log(pop)+splag| vote | (bzineq~whe) | BEZIRK, data = pooled))

# First stage ruggedness

summary(m91 <- felm(bzineq ~ rug| 0 | 0 | BEZIRK, data = pooled))

summary(m92 <- felm(bzineq ~ cath+ger+urban+log(pop)+rug+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage ruggedness

summary(m10 <- felm(yessh ~ 1| vote | (bzineq~rug) | BEZIRK, data = pooled))

summary(m11 <- felm(yessh ~ cath+ger+urban+log(pop)+splag| vote | (bzineq~rug) | BEZIRK, data = pooled))


texreg(list(m1,m2,m31,m32,m4,m5,m61,m62,m7,m8,m91,m92,m10,m11), 
       booktabs = T, dcolumn = T, use.packages = F, 
       no.margin = T, custom.header = list("OLS" = 1:2,"IV Soil" = 3:6,"IV Wheat" = 7:10, "IV Ruggedness" = 11:14),
       custom.model.names = c(rep("OLS",2), rep(rep(c("1. Stage","1. Stage","2. Stage","2. Stage"),1),3)),
       caption = "Landholding Inequality and Support for Suffrage Extension", label = "mainreg",
       custom.coef.map =  list("nut" = "Nutrient storage capacity", "whe" = "Climate Suitability Wheat","rug" = "Ruggedness",
                               "bzineq" = "Landholding Inequality","`bzineq(fit)`" = "Landholding Inequality","imr" = "Infant Mortality Rate",
                               "urban" = "Urbanization",  "log(pop)" = "ln Population",
                               "cath" = "Catholic District","ger" = "German District","splag" = "Spatial Lag","spbzineq" = "Spatial Lag", "(Intercept)" = "Intercept"),
       custom.gof.rows = list("Vote FEs" = c("Yes","Yes", rep(c("No","No","Yes","Yes"),3))),
       scalebox = 0.8, include.fstatistic = T, include.proj.stats = F, sideways = T , stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")


##### Table 2: Landholding Inequality and Support for Suffrage Extension: Social Control Mechanisms

m1 <- lm(sp_p_pup ~  spspendt + bzineq+cath+ger+urban+log(pop), bezsmall2)
m1se <- coeftest(m1,vcovHC(m1))[,2]
m1p <- coeftest(m1,vcovHC(m1))[,4]
m2 <- lm(teach_edu ~  sptedu + bzineq+cath+ger+urban+log(pop), bezsmall2)
m2se <- coeftest(m2,vcovHC(m2))[,2]
m2p <- coeftest(m2,vcovHC(m2))[,4]
m3 <- lm(edumean ~ spedu + bzineq+ cath+ger+urban+log(pop), bezsmall2)
m3se <- coeftest(m3,vcovHC(m3))[,2]
m3p <- coeftest(m3,vcovHC(m3))[,4]
m4 <- lm(votdens ~ spdens + bzineq+ cath+ger+urban+log(pop), bezall)
m4se <- coeftest(m4,vcovHC(m4))[,2]
m4p <- coeftest(m4,vcovHC(m4))[,4]
m5 <- lm(sect ~ spsect + bzineq+ cath+ger+urban+log(pop), bezall)
m5se <- coeftest(m5,vcovHC(m5))[,2]
m5p <- coeftest(m5,vcovHC(m5))[,4]


texreg(list(m1,m2,m3,m4,m5),  override.se = list(m1se,m2se,m3se,m4se,m5se),
       override.pvalues = list(m1p,m2p,m3p,m4p,m5p),
       booktabs = T, dcolumn = T, use.packages = F, 
       caption = "Landholding Inequality and Support for Suffrage Extension:  Social Control Mechanisms", label = "mainreg2",
       custom.coef.map =  list("bzineq" = "Landholding Inequality", "spspendt" = "Spatial Lag","sptedu" = "Spatial Lag",
                               "spedu" = "Spatial Lag","spdens" = "Spatial Lag","spsect" = "Spatial Lag",
                               "urban" = "Urbanization", "log(pop)" = "ln Population",
                               "cath" = "Catholic District","ger" = "German District","(Intercept)" = "Intercept"),
       custom.model.names = c("Education Spending", "Education Teacher","Education Performance","Voter Density", "Associations"),
       sideways = T, custom.note = "Heteroscedasticity-robust standard errors in parenthesis. %stars.", stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")

##### Table 3: Landholding Inequality and Support for Suffrage Extension: Alternative Mechanisms

m1 <- lm(turn_75 ~ sptu_75 + bzineq+ cath+ger+urban+log(pop), bezall)
m1se <- coeftest(m1,vcovHC(m1))[,2]
m1p <- coeftest(m1,vcovHC(m1))[,4]
m2 <- lm(turn_77 ~ sptu_77 + bzineq+ cath+ger+urban+log(pop), bezall)
m2se <- coeftest(m2,vcovHC(m2))[,2]
m2p <- coeftest(m2,vcovHC(m2))[,4]
m3 <- lm(excl ~  spexcl + bzineq+ cath+ger+urban+log(pop), bezsmall3)
m3se <- coeftest(m3,vcovHC(m3))[,2]
m3p <- coeftest(m3,vcovHC(m3))[,4]
m4 <- lm(outcan ~  spoutcan + bzineq+ cath+ger+urban+log(pop), bezsmall2)
m4se <- coeftest(m4,vcovHC(m4))[,2]
m4p <- coeftest(m4,vcovHC(m4))[,4]
m5 <- lm(imr ~ spimr+bzineq+ cath+ger+urban+log(pop), bezall)
m5se <- coeftest(m5,vcovHC(m5))[,2]
m5p <- coeftest(m5,vcovHC(m5))[,4]

texreg(list(m1,m2,m3,m4,m5),  override.se = list(m1se,m2se,m3se,m4se,m5se),
       override.pvalues = list(m1p,m2p,m3p,m4p,m5p),
       booktabs = T, dcolumn = T, use.packages = F, 
       caption = "Landholding Inequality and Support for Suffrage Extension: Alternative Mechanisms", label = "mainreg3",
       custom.coef.map =  list("bzineq" = "Landholding Inequality", "sptu_75" = "Spatial Lag","sptu_77" = "Spatial Lag",
                               "spexcl" = "Spatial Lag","spoutcan" = "Spatial Lag","spimr" = "Spatial Lag",
                               "urban" = "Urbanization", "log(pop)" = "ln Population","cath" = "Catholic District","ger" = "German District",
                               "(Intercept)" = "Intercept"),
       custom.model.names = c("Turnout 1875","Turnout 1877", "Excluded", "Non-Cantonal","Poverty"),
       sideways = T, custom.note = "Heteroscedasticity-robust standard errors in parenthesis. %stars.", stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")


############# Appendix

# Summary Statistics: Main Estimation

sumdat <- pooled %>% dplyr::select(whe,nut,rug,bzineq,yessh,cath,ger,excl,urban,pop) %>% st_drop_geometry() %>% mutate(pop = log(pop))

Nobs <- function(x) length(x)-sum(is.na(x))

datasummary(Heading("Yes Share")*yessh+Heading("Climate Suitability Wheat")*whe+Heading("Nutrient storage capacity")*nut+
              Heading("Ruggedness")*rug+Heading("Landholding Inequality")*bzineq + 
              Heading("Catholic District")*cath + Heading("German District")*ger + Heading("Share Disenfranchised")*excl + 
              Heading("Urbanization")*urban + Heading("ln Population Size")*pop ~ 
              Nobs + (Mean+Median+SD+Min+Max)*Arguments(na.rm=TRUE,fmt = "%.2f"),sumdat, fmt = "%.2f",
            title = "Summary Statistics: Main Estimation \\label{summain}")

# Summary Statistics: Mechanism

sumdat <- bezall %>% dplyr::select(KT,turn_75,turn_77,votdens,edumean,outcan,sect,sp_p_pup,teach_edu,cath,ger,excl,urban,pop,imr,bzineq)%>%
  st_drop_geometry() %>% mutate(pop = log(pop))


datasummary(Heading("Educarion Spending p. Pupil")*sp_p_pup + Heading("Share Poorly Trained Teacher")*teach_edu + 
              Heading("Education Performance")*edumean + Heading("ln Voter Density")*votdens + 
              Heading("Associations")*sect +Heading("Turnout 1875")*turn_75 +
              Heading("Turnout 1877")*turn_77 + Heading("Share Disenfranchised")*excl +Heading("Share Non-Cantonal Citizens")*outcan +
              Heading("Poverty (Infant Mortality)")*imr + Heading("Landholding Inequality")*bzineq + 
              Heading("Catholic District")*cath + Heading("German District")*ger + 
              Heading("Urbanization")*urban + Heading("ln Population Size")*pop ~ 
              Nobs + (Mean+Median+SD+Min+Max)*Arguments(na.rm=TRUE,fmt = "%.2f"),sumdat, fmt = 0,
            title = "Summary Statistics: Mechanisms \\label{summech}")


# Figure A5: Correlation Matrix of Instruments and Economic Development

bezall2 <- bezall %>% st_drop_geometry() %>% mutate(lpop = log(pop)) %>% rename("Wheat" = "whe", "Ruggedness" = "rug", "Soil" = "nut", "Urbanization" = "urban",
                                                                                "Infant Mortality" = "imr","Intra-Migration" = "intra")

ggpairs(bezall2[,c("Wheat","Ruggedness","Soil","Urbanization","Infant Mortality","Intra-Migration")],lower=list(continuous=wrap("smooth", method = MASS::rlm)), 
        upper = list(continuous = wrap("cor", size = 9))) +
  theme(strip.text = element_text(size = 30))

# Table A4: Non-linear specifications

summary(m1 <- felm(yessh ~ poly(bzineq,2,raw=T)| vote | 0 | BEZIRK, data = pooled))

summary(m2 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+poly(bzineq,2,raw=T)+splag| vote | 0 | BEZIRK, data = pooled))

# Second stage nutrient levels 

summary(m4 <-felm(yessh ~ 1| vote | (poly(bzineq,2,raw=T)~poly(nut,2,raw=T)) | BEZIRK, data = pooled))

summary(m5 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (poly(bzineq,2,raw=T)~poly(nut,2,raw=T))| BEZIRK, data = pooled))

# Second stage wheat

summary(m7 <- felm(yessh ~ 1| vote | (poly(bzineq,2,raw=T)~poly(whe,2,raw=T)) | BEZIRK, data = pooled))

summary(m8 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (poly(bzineq,2,raw=T)~poly(whe,2,raw=T)) | BEZIRK, data = pooled))

# Second stage ruggedness

summary(m10 <- felm(yessh ~ 1| vote | (poly(bzineq,2,raw=T)~poly(rug,2,raw=T)) | BEZIRK, data = pooled))

summary(m11 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (poly(bzineq,2,raw=T)~poly(rug,2,raw=T)) | BEZIRK, data = pooled))



texreg(list(m1,m2,m4,m5,m7,m8,m10,m11), 
       booktabs = T, dcolumn = T, use.packages = F, 
       no.margin = T, custom.header = list("OLS" = 1:2,"IV Soil" = 3:4,"IV Wheat" = 5:6, "IV Ruggedness" = 7:8),
       custom.model.names = c(rep("OLS",2), rep(rep(c("2. Stage","2. Stage"),1),3)),
       caption = "Landholding Inequality and Support for Suffrage Extension: Non-linear Specifications", label = "nonl",
       custom.coef.map =  list("nut" = "Nutrient storage capacity", "whe" = "Climate Suitability Wheat","rug" = "Ruggedness",
                               "bzineq" = "Landholding Inequality","poly(bzineq, 2, raw = T)1" = "Landholding Inequality", 
                               "poly(bzineq, 2, raw = T)2" = "Landholding Inequality Squared","`poly(bzineq, 2, raw = T).1(fit)`" = "Landholding Inequality", 
                               "`poly(bzineq, 2, raw = T).2(fit)`" = "Landholding Inequality Squared",
                               "imr" = "Infant Mortality Rate",
                               "urban" = "Urbanization", "excl" = "Share Disenfranchised", "log(pop)" = "ln Population",
                               "cath" = "Catholic District","ger" = "German District","splag" = "Spatial Lag","spbzineq" = "Spatial Lag", "(Intercept)" = "Intercept"),
       custom.gof.rows = list("Vote FEs" = c("Yes","Yes", rep(c("Yes","Yes"),3))),
       scalebox = 0.8, include.fstatistic = T, include.proj.stats = F, sideways = T , stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")

# Table A5: Determinants of Net Migration between Districts


summary(m1 <- lm(intra~ spintra +  bzineq + urban + imr + cath + +log(pop) + ger, bezsmall2))
m1se <- coeftest(m1,vcovHC(m1))[,2]
m1p <- coeftest(m1,vcovHC(m1))[,4]

texreg(list(m1), override.se = list(m1se),
       override.pvalues = list(m1p),
       booktabs = T, dcolumn = T, use.packages = F, 
       caption = "Determinants of Net Migration between Districts", label = "intramig",
       custom.coef.map =  list("spintra" = "Spatial Lag","bzineq" = "Landholding Inequality", 
                               "urban" = "Urbanization", "imr" = "Infant Mortality Rate","log(pop)" = "ln Population",
                               "cath" = "Catholic District","ger" = "German District",
                               "(Intercept)" = "Intercept"), 
       custom.note = "Heteroscedasticity-robust standard errors in parenthesis. %stars.", stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")

# Table A6: Proximate Determinants of Support for Suffrage Extension


summary(m1 <- lm(sufgem~ spsufgem +  sp_p_pup + teach_edu + edumean + votdens + sect + excl + outcan + imr, bezsmall))
m1se <- coeftest(m1,vcovHC(m1))[,2]
m1p <- coeftest(m1,vcovHC(m1))[,4]

summary(m2 <- lm(sufcan~ spsufcan + sp_p_pup + teach_edu + edumean + votdens + sect + excl + outcan + imr, bezsmall))
m2se <- coeftest(m2,vcovHC(m2))[,2]
m2p <- coeftest(m2,vcovHC(m2))[,4]

summary(m3 <- lm(yes_75~ spyes_75 + sp_p_pup + teach_edu + edumean + votdens + sect + excl + outcan + imr, bezall))
m3se <- coeftest(m3,vcovHC(m3))[,2]
m3p <- coeftest(m3,vcovHC(m3))[,4]

summary(m4 <- lm(yes_77~ spyes_77 + sp_p_pup + teach_edu + edumean + votdens + sect + excl + outcan + imr, bezall))
m4se <- coeftest(m4,vcovHC(m4))[,2]
m4p <- coeftest(m4,vcovHC(m4))[,4]

texreg(list(m1,m2,m3,m4), override.se = list(m1se,m2se,m3se,m4se),
       override.pvalues = list(m1p,m2p,m3p,m4p),
       booktabs = T, dcolumn = T, use.packages = F, 
       caption = "Proximate Determinants of Support for Suffrage Extension", label = "mech2",
       custom.coef.map =  list("spsufgem" = "Spatial Lag","spsufcan" = "Spatial Lag",
                               "spyes_75" = "Spatial Lag","spyes_77" = "Spatial Lag","sp_p_pup" = "Education Spending",
                               "teach_edu" = "Poorly trained Teachers", "edumean" = "Weak Educational Performance",
                               "votdens" = "Voter Density","sect" = "Associations","excl" = "Excluded","outcan" = "Non-Cantonal", "imr" = "Poverty",
                               "(Intercept)" = "Intercept"),
       custom.model.names = c("Vote I 1866","Vote II 1866", "Vote 1875", "Vote 1877"),
       sideways = T, custom.note = "Heteroscedasticity-robust standard errors in parenthesis. %stars.", stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")

# Table A7: Landholding Inequality and Support for Suffrage Extension (Disenfranchised as Control)


# OLS First Stage

summary(m1 <- felm(yessh ~ bzineq| vote | 0 | BEZIRK, data = pooled))

summary(m2 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+bzineq+splag| vote | 0 | BEZIRK, data = pooled))

# First stage nutrient levels  

summary(m31 <- felm(bzineq ~ nut| 0 | 0 | BEZIRK, data = pooled))

summary(m32 <- felm(bzineq ~ cath+ger+urban+excl+log(pop)+nut+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage nutrient levels 

summary(m4 <-felm(yessh ~ 1| vote | (bzineq~nut) | BEZIRK, data = pooled))

summary(m5 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (bzineq~nut) | BEZIRK, data = pooled))

# First stage wheat

summary(m61 <- felm(bzineq ~ whe| 0 | 0 | BEZIRK, data = pooled))

summary(m62 <- felm(bzineq ~ cath+ger+urban+excl+log(pop)+whe+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage wheat

summary(m7 <- felm(yessh ~ 1| vote | (bzineq~whe) | BEZIRK, data = pooled))

summary(m8 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (bzineq~whe) | BEZIRK, data = pooled))

# First stage ruggedness

summary(m91 <- felm(bzineq ~ rug| 0 | 0 | BEZIRK, data = pooled))

summary(m92 <- felm(bzineq ~ cath+ger+urban+excl+log(pop)+rug+spbzineq| 0 | 0 | BEZIRK, data = pooled))

# Second stage ruggedness

summary(m10 <- felm(yessh ~ 1| vote | (bzineq~rug) | BEZIRK, data = pooled))

summary(m11 <- felm(yessh ~ cath+ger+excl+urban+log(pop)+splag| vote | (bzineq~rug) | BEZIRK, data = pooled))

texreg(list(m1,m2,m31,m32,m4,m5,m61,m62,m7,m8,m91,m92,m10,m11),
       booktabs = T, dcolumn = T, use.packages = F, 
       no.margin = T, custom.header = list("OLS" = 1:2,"IV Soil" = 3:6,"IV Wheat" = 7:10, "IV Ruggedness" = 11:14),
       custom.model.names = c(rep("OLS",2), rep(rep(c("1. Stage","1. Stage","2. Stage","2. Stage"),1),3)),
       caption = "Landholding Inequality and Support for Suffrage Extension", label = "mainregwoexcl",
       custom.coef.map =  list("nut" = "Nutrient storage capacity", "whe" = "Climate Suitability Wheat","rug" = "Ruggedness",
                               "bzineq" = "Landholding Inequality","`bzineq(fit)`" = "Landholding Inequality","imr" = "Infant Mortality Rate",
                               "urban" = "Urbanization", "excl" = "Share Disenfranchised", "log(pop)" = "ln Population",
                               "cath" = "Catholic District","ger" = "German District","splag" = "Spatial Lag","spbzineq" = "Spatial Lag", "(Intercept)" = "Intercept"),
       custom.gof.rows = list("Vote FEs" = c("Yes","Yes", rep(c("No","No","Yes","Yes"),3))),
       scalebox = 0.8, include.fstatistic = T, include.proj.stats = F, sideways = T , stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\dagger")

# Figure A6: IV Estimates for the Individual Votes

ineqest <- matrix(NA, nrow = 12,2)

summary(m <- felm(sufgem ~ cath+ger+urban+log(pop)+spsufgem| 0 | (bzineq~nut) | 0, data = bezsmall))
ineqest[1,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[1,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sufgem ~ cath+ger+urban+log(pop)+spsufgem| 0 | (bzineq~whe) | 0, data = bezsmall))
ineqest[2,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[2,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sufgem ~ cath+ger+urban+log(pop)+spsufgem| 0 | (bzineq~rug) | 0, data = bezsmall))
ineqest[3,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[3,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sufcan ~ cath+ger+urban+log(pop)+spsufcan| 0 | (bzineq~nut) | 0, data = bezsmall))
ineqest[4,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[4,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sufcan ~ cath+ger+urban+log(pop)+spsufcan| 0 | (bzineq~whe) | 0, data = bezsmall))
ineqest[5,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[5,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sufcan ~ cath+ger+urban+log(pop)+spsufcan| 0 | (bzineq~rug) | 0, data = bezsmall))
ineqest[6,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[6,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_75 ~ cath+ger+urban+log(pop)+spyes_75| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[7,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[7,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_75 ~ cath+ger+urban+log(pop)+spyes_75| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[8,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[8,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_75 ~ cath+ger+urban+log(pop)+spyes_75| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[9,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[9,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_77 ~ cath+ger+urban+log(pop)+spyes_77| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[10,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[10,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_77 ~ cath+ger+urban+log(pop)+spyes_77| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[11,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[11,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(yes_77 ~ cath+ger+urban+log(pop)+spyes_77| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[12,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[12,2] <- m$se["`bzineq(fit)`"]

depvar <- c(rep(c("Vote 1866 I"),3), rep(c("Vote 1866 II"), 3), rep(c("Vote 1875"), 3), rep(c("Vote 1877"), 3))

instr <- rep(c("Soil","Wheat","Ruggedness"),4)

ivests <- data.frame(est = ineqest[,1], se = ineqest[,2], depvar = depvar, Instrument = instr)

ivests <- ivests %>% 
  mutate(depvar = factor(depvar, levels = c("Vote 1866 I","Vote 1866 II","Vote 1875","Vote 1877")))

ggplot(ivests, aes(x = depvar, y = est, colour = Instrument)) + 
  geom_point(position = position_dodge(width=.75)) + 
  geom_linerange(aes(ymin = est-se*1.645, ymax = est+se*1.645),position = position_dodge(width=.75)) + 
  geom_hline(yintercept = 0) + theme_bw()+ labs(x = "Dependent Variable", y = "Estimates with 90% CIs") +theme(axis.text=element_text(size=20),
                                                                                                               axis.title=element_text(size=26,face="bold"))

# Figure A7: Power Analysis: Sample Size Calculations


pool <- pooled %>% st_drop_geometry() %>% as.data.frame()

betas <- seq(-0.84, -0.36, 0.01)

samp <- matrix(NA, nrow = length(betas), ncol = 3)

for (i in 1:length(betas)){
  
  nut <- ivmodelFormula(yessh ~ cath + bzineq + splag + vote + urban + excl + ger + log(pop)|
                          cath + nut + splag + vote + urban + excl + ger + log(pop), data=pool)
  
  samp[i,1] <- IVsize(nut, power=0.8, type = "TSLS", beta = betas[i])
  
  whe <- ivmodelFormula(yessh ~ cath + bzineq + splag + vote + urban + excl + ger + log(pop)|
                          cath + whe + splag + vote + urban + excl + ger + log(pop), data=pool)
  
  samp[i,2] <- IVsize(whe, power=0.8, type = "TSLS", beta = betas[i])
  
  rug <- ivmodelFormula(yessh ~ cath + bzineq + splag + vote + urban + excl + ger + log(pop)|
                          cath + rug + splag + vote + urban + excl + ger + log(pop), data=pool)
  
  samp[i,3] <- IVsize(rug, power=0.8, type = "TSLS", beta = betas[i])
  
}


sampsiz <- data.frame(beta = betas, samp)

sampsiz <- sampsiz %>% rename(Soil = X1, Wheat = X2, Ruggedness = X3) %>% pivot_longer(!beta, names_to = "Instruments", values_to = "Sample")


ggplot(sampsiz, aes(y = Sample, x = beta, linetype = Instruments)) + geom_line() + theme_bw() + geom_hline(yintercept = 0 ) +
  labs(x = "Effect Size", y = "Necessary Sample Size") + 
  scale_y_continuous(breaks = seq(0, 2000, by = 200)) + theme(axis.text=element_text(size=20),
                                                              axis.title=element_text(size=26,face="bold"),
                                                              legend.title = element_text(size=20),
                                                              legend.text = element_text(size=20))


# Figure A8: IV Estimates for the Mechanisms  

bezsmall2 <- bezsmall2 %>% mutate(sp_p_pup2 = sp_p_pup/100)
bezall <- bezall %>% mutate(votdens2 = votdens/100, imr2 = imr/100)

ineqest <- matrix(NA, nrow = 30,2)

summary(m <- felm(sp_p_pup2 ~ spspendt + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezsmall2))
ineqest[1,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[1,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(teach_edu ~ sptedu + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezsmall2))
ineqest[2,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[2,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(edumean ~ spedu + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezsmall2))
ineqest[3,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[3,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(votdens2 ~ spdens + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[4,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[4,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sect ~ spsect + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[5,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[5,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sp_p_pup2 ~ spspendt + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezsmall2))
ineqest[6,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[6,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(teach_edu ~ sptedu + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezsmall2))
ineqest[7,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[7,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(edumean ~ spedu + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezsmall2))
ineqest[8,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[8,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(votdens2 ~ spdens + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[9,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[9,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sect ~ spsect + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[10,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[10,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sp_p_pup2 ~ spspendt + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezsmall2))
ineqest[11,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[11,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(teach_edu ~ sptedu + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezsmall2))
ineqest[12,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[12,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(edumean ~ spedu + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezsmall2))
ineqest[13,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[13,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(votdens2 ~ spdens + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[14,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[14,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(sect ~ spsect + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[15,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[15,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_75 ~ sptu_75 + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[16,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[16,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_77 ~ sptu_77 + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[17,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[17,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(excl ~ spexcl + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezsmall3))
ineqest[18,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[18,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(outcan ~ spoutcan + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezsmall2))
ineqest[19,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[19,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(imr2 ~ spimr + cath+ger+urban+log(pop)| 0 | (bzineq~nut) | 0, data = bezall))
ineqest[20,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[20,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_75 ~ sptu_75 + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[21,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[21,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_77 ~ sptu_77 + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[22,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[22,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(excl ~ spexcl + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezsmall3))
ineqest[23,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[23,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(outcan ~ spoutcan + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezsmall2))
ineqest[24,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[24,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(imr2 ~ spimr + cath+ger+urban+log(pop)| 0 | (bzineq~whe) | 0, data = bezall))
ineqest[25,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[25,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_75 ~ sptu_75 + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[26,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[26,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(turn_77 ~ sptu_77 + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[27,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[27,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(excl ~ spexcl + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezsmall3))
ineqest[28,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[28,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(outcan ~ spoutcan + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezsmall2))
ineqest[29,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[29,2] <- m$se["`bzineq(fit)`"]

summary(m <- felm(imr2 ~ spimr + cath+ger+urban+log(pop)| 0 | (bzineq~rug) | 0, data = bezall))
ineqest[30,1] <- m$coefficients["`bzineq(fit)`",1]
ineqest[30,2] <- m$se["`bzineq(fit)`"]

depvar <- c(rep(c("Edu_Spending","Edu_Teacher","Edu_Performance","Voter Density","Associations"),3), rep(c("Turnout 1875","Turnout 1877","Excluded","Non-Cantonal","Poverty"), 3))

instr <- rep(c(rep("Soil",5), rep("Wheat",5),rep("Ruggedness",5)),2)

ivests <- data.frame(est = ineqest[,1], se = ineqest[,2], depvar = depvar, Instrument = instr)

ivests <- ivests %>% 
  mutate(depvar = factor(depvar, levels = c("Edu_Spending","Edu_Teacher","Edu_Performance","Voter Density","Associations","Turnout 1875","Turnout 1877","Excluded","Non-Cantonal","Poverty")))

ggplot(ivests, aes(x = depvar, y = est, colour = Instrument)) + 
  geom_point(position = position_dodge(width=.75)) + 
  geom_linerange(aes(ymin = est-se*1.96, ymax = est+se*1.96),position = position_dodge(width=.75)) + 
  geom_hline(yintercept = 0) + theme_bw()+ labs(x = "Dependent Variable", y = "Estimates with 95% CIs") +theme(axis.text=element_text(size=20),
                                                                                                               axis.title=element_text(size=26,face="bold"))